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We present a new derivation of the perturbation equations governing the oscillations of relativistic 
non-rotating neutron star models using the ADM-formalism. This formulation has the advantage 
that it immediately yields the evolution equations in a hyperbolic form, which is not the case for the 
Einstein field equations in their original form. We show that the perturbation equations can always 
be written in terms of spacetime variables only, regardless of any particular gauge. We demonstrate 
how to obtain the Regge- Wheeler gauge, by choosing appropriate lapse and shift. In addition, not 
only the 3-metric but also the extrinsic curvature of the initial slice have to satisfy certain conditions 
in order to preserve the Regge- Wheeler gauge throughout the evolution. We discuss various forms of 
the equations and show their relation to the formulation of Allen et al. New results are presented for 
polytropic equations of state. An interesting phenomenon occurs in very compact stars, where the 
first ring-down phase in the wave signal corresponds to the first quasinormal mode of an equal mass 
black hole, rather than to one of the proper quasinormal modes of the stellar model. A somewhat 
heuristic explanation to account for this phenomenon is given. For realistic equations of state, the 
numerical evolutions exhibit an instability, which does not occur for polytropic equations of state. 
We show that this instability is related to the behavior of the sound speed at the neutron drip point. 
As a remedy, we devise a transformation of the radial coordinate r inside the star, which removes 
this instability and yields stable evolutions for any chosen numerical resolution. 



I. INTRODUCTION 



The driving force behind the theoretical studies in gravitational radiation is to predict and investigate the sources 
of gravitational waves that will be observed by interferometric (GEO600, LIGO, VIRGO, TAMA, LISA) and bar 
(ALLEGRO, AURIGA, EXPLORER, NAUTILUS) detectors. The goal of these studies is not limited to facilitating 
the detection of signals but also to provide tools to extract as much information as possible about the physical nature 
of the sources. Neutron stars have been pointed out as obvious source candidates, since one channel by which an 
oscillating neutron star looses energy is via the emission of gravitational radiation. 

This radiation basically consists of a superposition of its characteristic oscillatory modes, which can be grouped into 
two families: i) The fluid modes jl]|2), which have a Newtonian counterpart, but which through their coupling to the 
spacetime are damped because they now can emit gravitational waves, ii) The spacetime modes which have no 

Newtonian counterpart, and which couple only weakly to the fluid. (In the odd parity case, they do not couple to the 
fluid at all). These modes usually are strongly damped, however, for ultra-relativistic stars, the spacetime curvature 
can be so strong that it can trap impinging gravitational waves. Those "trapped" modes HQ are quite long lived, 
since they only leak out slowly from inside the gravitational potential well created by the neutron star. For recent 
comprehensive reviews on oscillation modes of neutron stars and black holes, see ^,|). 

However, just knowing the different oscillation modes of a neutron star is not enough. To be able to detect and to 
interpret a signal in a gravitational wave detector, it is crucial to have accurate templates of the waveforms that result 
from astrophysical events, such as the formation of the neutron star after the collapse of a progenitor star at the end 
point of nuclear burning. To that purpose, two ingredients are necessary. Firstly, one has to have initial data that 
represent a realistic stage of the astrophysical system under consideration. Secondly, one needs a numerical evolution 
code that can evolve these initial data in a stable and accurate way. 

It is only by fully nonlinear simulations that accurate waveforms of oscillating neutron stars after core collapse or 
binary merger can be obtained. Therefore it is a major goal of several groups to push the further development of 
nonlinear evolution codes. The biggest obstacle, however, that still prevents one from obtaining accurate results, is 
the enormous computational expenditure to solve the full set of Einstein equations in 3D. Even modest resolutions 
easily exceed the capabilities of the largest and fastest present day computers. 

It is therefore crucial to have some alternative methods which require less computing power but still give valuable 
physical insight. One of them is perturbation theory. In the non-rotating case, due to the spherical background, 
one can completely separate off the angular dependence, and the evolution problem boils down to the integration of 
ID wave equations. Slowly rotating stars can still be treated in ID as long as the deformation due to rotation is 
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negligible. But for rapidly rotating neutron stars, the spherical symmetry is broken, and one can only separate off 
the azimuthal angular dependence. Hence the problem becomes 2D, which is nevertheless still numerically tractable. 
The results from those linear evolutions are expected to give accurate waveforms in the range where the neutron star 
oscillates only weakly. Furthermore, they can be used as testbeds for the nonlinear calculations. Indeed, recent results 
of nonlinear evolutions have confirmed the existence of the various modes predicted by linear perturbation theory p0[ . 

The evolution of polar perturbations of non-rotating stars has been first addressed by Allen et al. jy]. In their 
formulation they describe the oscillations with three coupled wave equations, two for the metric perturbations in and 
outside the star, and one for the fluid inside the star. The evolution of various initial data have shown that generic 
initial data can indeed excite some of the stellar modes such as the /-mode, some p-modes and the first w-mode. 

In this paper we rederive the perturbation equations using the (3+1) or ADM- form [JL2| of the Einstein equations. 
These equations are much better suited for the numerical initial value problem, since they immediately yield a 
hyperbolic evolution system provided that a suitable gauge has been chosen, whereas the original field equations 
contain mixed time and space derivatives. It is only by the introduction of new variables that those equations can 
be cast into a well suited form. The ADM- formalism has already been addressed by Moncrief jl3],[l4]], where it was 
used to derive a Hamiltonian gauge invariant formulation of the perturbation equations. Albeit preferable from a 
conceptual point of view, since they are independent of any gauge, those equations do not prove particularly useful for 
numerical evolutions. Our starting point therefore are the field equations written as a set of evolution equations for 
the metric and the extrinsic curvature of a 3-dimensional spatial hypersurface, together with the constraint equations, 
which have to be satisfied at every instance of time. 

As an immediate consequence of using the ADM-formalism, it follows that it is always possible to eliminate the fluid 
variables and thus to write down the evolution equations in terms of metric quantities only. This has already been 
shown by Chandrasekhar and Ferrari for the diagonal gauge and by Ipser and Price [|l6| for the Regge- Wheeler 
gauge. However, from our formalism it is clear that this feature is independent of the chosen gauge. 

After having expanded the equations in spherical harmonics, we will choose the Regge- Wheeler gauge as the basis 
of further investigations. Due to the spherical harmonics the equations are independent of the azimuthal order m and 
can be divided into two uncoupled sets according to their behavior under parity transformation. The polar or even 
parity equations transform as (—1)', whereas the axial or odd parity equations as (— Our main focus will be 
the polar perturbations, nevertheless, due to their simplicity, we also present the axial equations in the ADM-form. 
The derivation of the polar equations is much more involved, because the "raw" forms lead to numerical instabilities 
at the origin due to indefinite expressions, which is a well known consequence of using spherical coordinates. We are 
therefore forced to recast the polar equations into a form that is well behaved at the origin. 

Having forged the equations into a form that yields stable evolutions for polytropic equations of state, we face a 
new instability close to the stellar surface when we try to switch to realistic equations of state. In this paper we use 
a quite recent equation of state, called MPA fr^| . Yet, this instability is not a consequence of just this particular 
equation of state, it will occur for any realistic equation of state. As we shall show, the instability results from the 
behavior of the sound speed at the neutron drip point, where the equation of state suddenly becomes very soft. This 
is accompanied by a sharp drop of the sound speed, which can give rise to a numerical instability if the spatial grid 
resolution is chosen too low. However, by increasing the resolution this instability will eventually vanish. In the non- 
rotating case, having to resort to higher resolutions does not present too great an obstacle, since the equations are 
ID wave equations. However, this instability presumably also occurs for rotating stars, where the equations are 2D, 
and the required resolution in order to obtain stable evolutions might stretch the computing times to unacceptably 
large values. Therefore we devise a coordinate transformation of the radial coordinate inside the star, which enables 
stable evolutions for any spatial resolution. 

The remainder of this paper is organized as follows: In Sec. II we derive the relevant equations. Section III is 
devoted to a discussion of the boundary conditions. In Sec. IV we briefly comment on the construction of initial 
data, and in Sec. V we present results for polytropic equations of state. Finally, in Sec. VI we discuss the instability 
associated with the use of realistic equations of state, and we present the coordinate transformation that removes it. 
Conclusions are briefly presented in Sec. VII. 

We adopt the metric signature ( — | — I — (-) and work in geometric units with G = c = 1. Spacetime and spatial 
indices are denoted by Greek and Latin letters, respectively. Derivatives with respect to the radial coordinate r are 
sometimes denoted by a prime and derivatives with respect to the time coordinate t by an overdot. 

II. DERIVATION OF THE PERTURBATION EQUATIONS 
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A. The unperturbed stellar model 



The background geometry of a non-rotating, spherically symmetric star is given by the line element 

ds 2 {0) = -e 2l/ dt 2 + e 2X dr 2 + r 2 {d0 2 + sin 2 6d(f> 2 ) , (1) 

where v and A are functions of the radial coordinate r. We model the star as a perfect fluid whose energy-momentum 
tensor has the form 

Tfiu = (p + e)u fi u v +pg^i, , (2) 

with p denoting the pressure, e the energy density, and the 4-velocity of the fluid. In the fluid rest frame 
has only one non- vanishing component u° — e~ u . Einstein's equations G M „ = 87rT M „ and the conservation equations 
T^ v . v = yield the following three independent structure equations for the four unknown A, v, p and e: 

i _ P 2A 

A' = + 47rre 2A e (3a) 

2r 
e 2A _ 1 

v 1 = — V 4irre 2X p (3b) 

p' = -v'(p + e). (3c) 

We will refer to those equations as Tolman-Oppenheimer-Volkoff (TOV) equations, although they are usually written 
in terms of the gravitational mass function m, which is related to A through 

o\ 2m . , 

e- 2X = 1 . (4) 

r 

To fully determine this system of equations, an equation of state must be supplemented. In this paper, we restrict 
ourselves to barotropic equations of state, where the pressure is a function of the energy density alone: 

p = p(e). (5) 

In particular, we will use either a polytropic equation of state in form of 

p = ne r (6) 

or a tabulated realistic zero temperature equation of state. To obtain the stellar model, we have to integrate the TOV 
equations (^|) together with the equation of state (^|) from the center up to the point where the pressure p vanishes. 
This then defines the surface R of the star. From Birkhoff theorem it follows that the exterior vacuum region of the 
star is described by the Schwarzschild metric with the mass parameter M = m(R). 



B. The perturbation equations 

1. The general form 

Let us recall that the general line clement, when written in the ADM- form, reads 

ds 2 = - (a 2 - [3 k [3 k ) dt 2 + 2$ dt dx i + 7^ dx l dx j , (7) 

where a, (3 k , and 7y are lapse function, shift vector, and spatial metric, respectively. Comparing the ADM-metric (Q) 
with the background metric ([!]) reveals that the background shift vanishes and the background lapse function is 
a (o) = e u ■ In addition, the extrinsic curvature of constant t slices vanishes, since the background metric is static. 

Let a and f3i now be the perturbations of the lapse function and the (covariant) shift vector, and let hij denote the 
perturbations of the spatial metric. The line element describing the perturbations to the background metric (Q) then 
reads 

ds 2 (1) = ~2e"a dt 2 + 2/3, dt dx l + hy dx l dx j . (8) 
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The perturbations of the energy-momentum tensor (Q) can be written in terms of 8e, 8p, and Su^, which are the 
(Eulerian) perturbations of the energy density e, pressure p, and (covariant) 4- velocity it M , respectively. We assume 
the equation of state (||) to hold for both the unperturbed and perturbed configurations (isentropic perturbations), 
therefore the perturbations 5e and Sp are related through 

§p= ?P Se = lse = C*6e, (9) 
ae e 

where C s is the sound speed inside the fluid. 

The dynamic equations that govern the metric and extrinsic curvature perturbations hij and fey read after lin- 
earization 



dthij = (3 k d kllJ + lkl d 3 p k + lkj d^ k - 2e lJ k lJ (10) 

(11) 



dtkj = -didja + T%d k a + ST%d k e" + a + 4tt (p - e) 7y ] 



+ e v [ SRij + 4tt (0 - e) hy + Se (C 2 S - l) Ty ) 



Therein, the 8T k ^ represent the perturbations of the spatial Christoffel symbols T\j 

ST % = (hu-,j + h w - h ij; i) 

I (12) 
- ^l kl (djhu + dihtj - dthii) ~ h lml kl T™ 

and SR^ the perturbations of the spatial Ricci tensor Rij 

6R ij =6T% ;k -6T k ik . tj 

= d k ST% d 3 5Y k lk - 6T k mk T™ + 5T™r k mk - 5T k m] T™ k - 5T™ k r k mj . 

We should note that the indices of perturbations are lowered or raised with the background 3- metric 7-y. In con- 
structing initial data, the perturbations must satisfy the linearized version of Einstein's constraint equations, namely 

j ij SRij - h ij Rij = KinSe (14) 
^ k (dik jk - d k lk - T l lk k ol + T l jk k u ) = -8n{p + e)5u l . (15) 

From the linearized energy-momentum conservation equations 5T^. U = 0, we can deduce equations of motion for the 
fluid perturbations 5e and Siii, which complete the set of dynamical equations. 

However, we can as well dispense with the fluid equations, since it is possible to eliminate the energy density 
perturbation Se from the evolution equation ( |TT| ) by virtue of the Hamiltonian constraint (|l4|). In this way we can 
obtain a consistent system of evolution equations for the metric and extrinsic curvature perturbations alone. The 
constraints then can serve as a means to compute the matter perturbations be and Su.i. We should stress that the 
possibility to completely describe the oscillations of neutron stars with spacetime variables only is a general feature 
of the perturbation equations and does not depend on any specific gauge choice. 

Due to the spherical symmetry of the background, we can eliminate the angular dependence of the perturbation 
equations by expanding them into spherical tensor harmonics. Those tensor harmonics can be divided into two subsets 
that behave differently under parity transformation. Under space reflection the even parity or polar harmonics change 
sign according to (—1)', whereas the odd parity or axial harmonics transform like (— 1) + . Here, I is the number that 
labels the multipole of the spherical harmonics Y/ m . Because of the spherical symmetry, the axial and polar equations 
are decoupled and degenerate with respect to the azimuthal number m. 

Without picking a specific gauge, the metric perturbations can be expanded using the tensor harmonics given first 
by Regge and Wheeler in the following way (symmetric components are denoted by an asterisk) : 

00 1 

a = Yl e "3f%» ( 16a ) 

1=2 m=-l 
00 I 

(3k = E (s l 2 m Ylrnrt m ^™ + ^ m $L m ) (16b) 



1=2 m=-l 



h — I I (If, \ 



1=2 m=-l 
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where 



$L m =(-sin- 1 ^, sin^y ;m , (18) 



and 



*!S= n ^™ ( 19 ) 



1 

^ ~ I sin 2 

*s=(£; < 2o » 

xS-"«(-*5L'*-!£:). pi) 



with 



lf| m = -h^- cot0— - sin"" 9— Yi 



2 \d6 2 86 8<f>- 

1 8 2 \ 



(22) 



- y im=(4-cotfl > )^ m . (23) 



89 ) d<p 

The notation has been chosen such that the coefficients Si represent the scalar parts of h^ u , namely a,f3 r , and h rr , 
whereas the Vi stand for the vector components fio,P<t» We and h r(j> . Lastly, the % represent the tensor components 
hgg^hetf, and h^. Note that this expansion includes both polar and axial harmonics. The latter are represented by 
<I>^™ and with respective coefficients V2, Va and T3. Similarly, the extrinsic curvature tensor kij will be expanded 



as 



(24) 

i= 2m =-l\ * r V 4 V «/3 +A 5 + X Q/3 J y 

Here, K l ™ and ifg m are the axial coefficients. 
Last not least, we need the matter variables 

00 / 

Se = E P lmY ^ ( 25 ) 

/=2 m=-/ 
00 I 

ui = E E ( fl ' my im, 4 m *L m + 4 m < m ) • (26) 

i=2 m=-; 

The sum over the multipoles starts from 1 = 2, since we are only interested in perturbations which are associated 
with the emission of gravitational radiation. Besides, the non-radiative multipoles I = and I = 1 would have to be 
treated in a different way. 

Using the above expansions, we obtain 12 evolution equations for the coefficients of the 3-metric hij and the extrinsic 
curvature kij. However, it is clear that this set cannot be used for the evolution, for we have not specified any gauge 
yet. Within perturbation theory, picking a specific gauge usually means to define a gauge vector jf and transform 
the metric perturbations according to 

KT ^ <t v** ^ • ( 2? ) 

In our case, however, we follow the spirit of the ADM- formalism, that is we pick the gauge by prescribing the 
coefficients of lapse a and shift Pi, for which we do not have any evolution equations. Yet, this is not enough as we 
shall see. To fully fix the gauge, we also have to impose certain constraints on the initial data. 
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The most common gauge is the Regge- Wheeler gauge |Jj|, which can be obtained by setting a certain number of 
metric perturbation coefficients to zero. Translating this into the ADM-formalism means to choose shift and lapse in 
such a way that for initial data which obey the Regge- Wheeler gauge, the evolution will preserve it. 

Regge and Wheeler have chosen their gauge such that the metric coefficients V x , V3 m , Ty , and T^ m vanish. By 
choosing shift and lapse as (from now on we omit the indices I and m) 

51 = -\S 3 (28) 

5 2 = 2e u K 2 (29) 
Vi=0 (30) 
V 2 = e»K 6 , (31) 

we obtain the following evolution equations for the coefficients Vi,Ti and T3: 

§-V 3 = (32) 

-|fi = -2e v K 4 (33) 

|f 3 = . (34) 

In addition the evolution equation for the extrinsic curvature component A4 depends only on T\ and V3. Thus, only 
for initial data satisfying V3 = T\ = T3 = and A4 — 0, the evolution equations guarantee the vanishing of those 
coefficients throughout the whole evolution. 

Hence, the Regge- Wheeler gauge does not only impose constraints on the initial metric perturbations but also on 
the extrinsic curvature perturbations. The vanishing of the coefficient A4 is a crucial feature of the Regge- Wheeler 
gauge and has to be imposed on the initial data. This is in contradiction with the statement of Andersson et al. ]l9| , 
where the authors claim that there are no constraints on the odd or even parity nature of the extrinsic curvature. This 
is not true, since the vanishing of A4 must be ensured on the initial slice, otherwise the evolution cannot be performed 
in the Regge- Wheeler gauge. For data with non- vanishing A4, we would have dT±/dt 7^ 0, which is incompatible with 
the Regge- Wheeler gauge, because it would lead to a non-vanishing T\ even if T\ was initially set to zero. Therefore 
the initial data presented in for two colliding neutron stars are not adequate for an evolution code which makes 
use of the Regge- Wheeler gauge, since the coefficient ki in Eqs. (82) of (which corresponds to our A4) does not 



vanish. We will discuss the issue of constructing initial data in more detail in section IV. 

Having switched to the Regge- Wheeler gauge, we are left with evolution equations for the three metric variables 
S3, V4, and T2, and all the extrinsic curvature variables Ki save K4. (In the more common notation of Regge and 
Wheeler it is S3 = H 2 , V4 = h\ and T 2 = A"). We now split the equations with respect to their behavior under 
parity transformation. Let us first focus our attention to the axial case. 



2. Axial perturbations 



If in the expansions (16a) and (|24|) we use new variables 



M — X-tr 

e 14 



A 3 = 2e- A A 3 
A 6 = r 2 e A A 6 , 



we obtain the following set of evolution equations: 



dV<t _ 2u-2\ 

dt 

dA 3 

dt 
dK 6 

dt 



1(1 + 1) 

dVj 
dr 



dK 6 
dr 

- 2 



/ - A' - - A, 



e 2A A 3 



• Va 



(35a) 
(35b) 
(35c) 

(36a) 

(36b) 
(36c) 
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and one constraint equation: 



From the conservation law 5T^ V . U = 0, it follows that 

£-0, (38) 

which means that axial gravitational waves do not couple to the stellar fluid. Of course, this does not mean that 
there cannot be any axial modes at all, in fact, there exists a whole spectrum of axial spacetime modes, which is quite 
similar the polar spectrum p0| . It is only the fluid modes that are missing. 

The three evolution equations (|3^) can be combined to yield a single wave equation for the metric perturbation V4, 
or better for the variable 

Va 

Q ■■= — , (39) 
r 

which completely describes the dynamical evolution of axial perturbations: 

d 2 Q d*Q 2v (, . , 6m lQ±m n r ,„. 

Herein, denotes the tortoise coordinate defined by 

^ = ^ ■ (41) 
dr* 

In the exterior, where p and e vanish and m — M, Eq. ( f40| ) reduces to the well known Regge- Wheeler equation 



3. Polar perturbations 

For the polar case we have only two dynamical metric perturbations S3 and T 2 , since the lapse S± is proportional 
to S3, and the only non- vanishing component S 2 of the shift is proportional to K 2 . To obtain the relevant equations, 
we again choose a somewhat different expansion. For the metric we use 

V (-+rS)Y lm (42) 



2 V r 

„2A 



A - (e 2X K 2 , 0, 0) Y lm (43) 
e 2A(Z +r5 ) 

hij = ( rT \Y lm , (44) 



Or sin 2 T 



and for the extrinsic curvature 



e 2A Xi -re 2X K 2 JL -re 2X K 2 JL 
h 3 = I * r 2 {K 5 - 2K 2 ) I Y im ■ (45) 

1 * r 2 sin 2 6> (Jf 5 - 2K 2 ) 



The expansion of the matter variables reads 

5e = P - Y lm (46) 
r 

^ ( u 2 d \ 

0Ui = \ui , u 2 — , u 2 — Y im . (47) 



r 



Of course, these expansions are to be understood as sums over all I and m, analogous to the expansions (16a) and 
(p||). It is only through this particular choice of expansion coefficients (e.g. in a, h rr and kgg) that we obtain a system 
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of equations that can be numerically integrated in a quite straightforward and - this is the main point - in a stable 
way. 

Another consequence of this decomposition is that it gives some physical meaning to the metric perturbations 5 
and T. From Eq. (0), we can see that T represents the conformal part of the spatial perturbation hij, for if we set 
5 = 0, the metric is conformally flat. Obviously, 5 itself then represents the deviation from conformal flatness. 

Finally, we should note that this decomposition is quite similar to the one used by Allen et al. Jll]]. Their variables 
F and SaiUu are related to ours as follows: F = T and Saiuu — e 2 "S. 

In order to avoid numerical problems at the origin, we have to replace K\ by the following quantity 



r l K := K x + 2r 



\ dr 



X'K, 



In this way, we obtain a system of five coupled evolution equations, which are of first order in time and second order in 
space. There are two equations for the metric variables 5 and T and three more for the extrinsic curvature variables 
K, K 2 and K 5 : 



dt 
dK 
~dt 



2v~2\ 



d 2 S 



OS 



+ (5l /-A') — + 4(^ + 5- 



dr 2 



dr 



3* 

r 



„2A 



- 1 



, 2A 



1(1+1) 



s 





(49a) 
(49b) 



dT 
~di 

dK 5 
dt 

dK 2 
dt 



„2v-2\ 



d 2 T 
dr 2 



dT 
or 



X' 



J2\ 



,,2A 



1(1+1) 



T + 2{rv' + rX' - 1) S 



2v-2\ 



+ %ire 2l/ (l-C 2 ) p 
OS 



1)5 + 2— T 



(49c) 
(49d) 

(49e) 



We could easily convert this system into a first order system both in time and space by addi ng an oth er tw o evolution 
equations for the first derivatives 5' and T'. However, the form of the first four equations (49a) - (49d), which are 
independent of K 2 , suggests to rather convert them into two coupled wave equations for 5 and T 



d 2 s 
dt 2 



d 2 r 
W 



„2i/-2A 



d 2 S 



dS 



dr 2 



dr 



X' 



„2A 



^ + (5^'-A')^+ 4« + 5-+3--2 



„2A 



1(1 + 1) 



5 



2u-2\ 




T + 2 (rv' + rX' -1)5 



(50) 



(51) 



(i - C 2 S ) p , 



which are equivalent to Eqs. (14) and (15) of Allen et al. |TT[. As can be seen, the wave equation for 5 ( |50|) is totally 
decoupled from the fluid variable p, which only couples to the metric perturbation T in Eq. (|5l]). The equation for 
K 2 is only necessary in the interior region, where it couples to the hydrodynamical equations, which follow from 
energy-momentum conservation ST 11 "^ = 0, and are given by 



dp 
dt 



dux 
~di 
du 2 
~dT 



a 2v-2X 

+ (P- 
^2°P 



diii 1 du 2 



r dr 

dK 2 

1- 

dr 

l + C 2 ) 



^--^+(3^-A'+-]ux-(3 



"2 



(2 + rX')K 2 - 7 ^K- i -K b 



re K 2 



(C 2 )' 



1 



dS dT 



+ ar +2rS 



C 2 sP -^(p + e) {r 2 S + T) 



(52a) 

(52b) 
(52c) 
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Here, we have defined Ui :— (p + e) Uj. By introducing the enthalpy perturbation 



H 



p + e 



-P 



(53) 



the fluid equations assume a more convenient form: 



dH 
~dt 



b 2v-2\q2 



dui 



—— + (2z/ - A') mi + 2--e 



dr 



r r 
3 



"2 



V \U\ u 2 

r 



dK r" 
+ C 2 ( r-^- + (2 + r\') K 2 -—K- ^K 5 ) - rv' ' K 2 



dm dH 1 / 2 dS dT 

= ^ ^ r- 2 — + — + 2rS 



dt dr 



dr dr 



du 2 



H--(r 2 S + T) 



(54a) 

(54b) 
(54c) 



Interestingly, from Eqs. ( |54b| ) and (54c), it follows that the coefficients u\ and u 2 are not independent of each other 
but rather are related via 



du 2 , , 



(55) 



where F is a time independent function, which has to be fixed by the initial data. The above system (|5J), too, can be 
cast into a second order wave equation for H, which is equivalent to Eq. (16) of Allen et al. |ll| (the different signs 
in the terms containing 5* and T are correct): 



d 2 H 



2v-2\ 



c\ 



d 2 H 
dr 2 



+ (C 2 (2,' - A') - „') ^ 4 ( C; ( — + 4 



A' 2 J(i + m , y , A' 



H 



I r r 



— 2rv' + - 



(r 2 S + T) 



(56) 



As already mentioned in subsection B.l, we do not necessarily need the fluid equations, for we can eliminate them by 
means of the Hamiltonian constraint, which relates the fluid variable p to the two metric variables S and T: 



8^p = -0 + A'f + + (2 - 2rA' + |e»>l(I + 1 , ) s 



,2A 



1 .A' 



„2A 



T . 



(57) 



With p substituted by the Hamiltonian constraint (|57j), the wave equation for T ( |5l| ) in the interior now reads 



d 2 T 

u 1 2u-2\ n 2 
2v-2\ 



d 2 T dS x ,dT 



A'— + 2rA' - 2 - -e 2X l(l + l))S + 
or \ 2 / 

V e 2A - 1 



<9r 2 <9r 



-,2A 



1 .A' 



,2A 



T 



T 



(58) 



In the exterior, both S and T propagate with the local speed of light e u , in the interior, however, T changes its 
character and propagates with the local speed of sound e v ~ x C s . 

The last set of equations that is still missing are the momentum constraints, which link the velocity perturbations 
to the extrinsic curvature variables: 



167re 2l/ (p + e) u\ 



- 2- 



dr dr 
dr 



tK- (3i/ + 3A'-e 2A ^-±^ 
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(59a) 
(59b) 



Those constraint equations have to be solved in order to obtain physically valid initial data. We postpone the 
discussion on how to construct initial to section IV . 







With p being eliminated, S and T in the interior now become independent variables, and the Hamiltonian constraint 
( |57j ) serves as a definition for p. In the exterior, however, S and T are not independent but have to satisfy the 
Hamiltonian constraint with p set to zero. Unfortunately, we cannot use the Hamiltonian constraint to further 
eliminate one of those variables, but it is possible to combine S and T to form a new variable Z (we use the definition 
(20) of Allen et al. @) 



2e 2 



(9 A/f 
— - 2-1(1 + 1) ) T - 2r 2 S 



(60) 



with 



A = Z(/ + l)-2 + — . (61) 



r 



Z then satisfies a single wave equation, the famous Zerilli equation that was first derived in 1970 by F. Zerilli pifl in 
the context of black hole oscillations: 

d 2 Z _ d 2 Z 2l/ n 2 (n + l)r 3 + 3n 2 Mr 2 + 9nM 2 r + 9M 3 

'W^W 2 ^ C " r 3 (nr + 3M) 2 ' ( ' 

Here, we use 2n = 1(1 + 1) — 2, and is again the tortoise coordinate defined in Eq. (|4l|). 
It is also possible to invert Eq. ( |60| ) and express S and T in terms of Z through 

T = re 2v Z' + ( h(l + 1) - ^-e 2 ^ Z (63a) 



2 y ' rA 



2 N 

e 2v 



(63b) 



Lastly, we should note that the radiated energy at infinity can be computed from p2] | 
dE 1 (j + 2)! 2 

^-^2^jr^y\ Z ^\ ■ (64) 

l,m y ' 

Since the Zerilli function contains the full polar gravitational wave information, it should be possible to recover the 
original metric and extrinsic curvature variables from it. For S and T we have already given the relevant formulas in 
Eqs. (|6^). Since K = S and K5 — T, they also can computed from (|6q ) with Z replaced by Z. In the exterior the 
momentum constrain ts ( p3) can be combined to give an algebraic relation for Finally, the last remaining quantity 
Ki follows from Eq. (|48|). 

In the derivations of all the above equations we have made extensive use of the computer algebra program MAPLE 
V in order to avoid possible mistakes. Furthermore, we have checked that our equations and our numerical results 
agree with those of Allen et al. |0|. 



III. BOUNDARY AND JUNCTION CONDITIONS 



There are three boundaries we have to take care of: The origin, the surface, and the outer boundary of the numerical 
grid, which lies somewhere outside the star. 

At the origin r = we have to demand all variables to be regular. From Taylor expansion around r — 0, we can 
infer the analytic behavior of the various variables. Close to the origin, S and T, for instance, are both proportional 
to r l+1 . 

At the outer boundary, we impose the Sommerfeld boundary condition, i.e. we require the waves to be purely 
outgoing. If computational time does not matter, we can even put this boundary so far away that any contamination 
that enters the grid from there has not enough time to travel to the region where we extract the signal. 

The third boundary is the surface of the star at r = R, which is formally defined by the vanishing of the total 
pressure P. Since the perturbations will slightly deform the star, the perturbed surface will be displaced by an amount 
£ l with respect to the unperturbed location at r — R. If the coordinates of the unperturbed surface are denoted by x % Rl 
the vanishing of the total pressure P at the displaced surface translates to P(t, x % R + £ J ) = 0. From Taylor expansion 
to first order we find that the (Eulerian) pressure perturbation Sp at the surface has to obey 

*P = -CP ■ (65) 
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Unfortunately, this is not a very convenient boundary condition, since we neither use Sp nor the displacement vector 
in our set of evolution equations. Therefore we must relate this condition to the variables we use. We will try 
to find a condition that gives us the time evolution of 5e at the stellar surface. The first step is to use the relation 
8e = p' / e' Sp, which gives us 



d_ 

m 



5c 



d 

-e'—e 

or 



The time derivative of £ r can then be related to the r-component of the 4-velocity u r [£3[ 
r)f r 

= e~ 2X (e»6u r - fa) . 
After expansion in spherical harmonics, we finally obtain 



r=R 



The equivalent equation for the quantity H as defined in (p3|) reads 



rK 2 + e 2l, -' 2X 



u 2 

ui 

r 



(66) 



(67) 



(68) 



(69) 



r=R 



Incidentally, this expression can be obtained directly from the evolution equation ( p4a[ ) just by setting C s to zero. In 
the same mann er, we can obtain the boundary condition for the wave equation (pq). Also, Eq. (p8|) follows immediately 



from Eq. ( 52a ), if one sets p = e = p' = 0. 



For poly trop ic equ ation s of state, it is always p = e = C s = p' = at the surface of the star, hence the evolution 
equations (52a) and (54a) automatically lead to the right boundary conditions. For realistic equations of state, the 
sound speed at the surface should be that of iron, which is very small compared to the sound speed inside the core, 
where it might reach almost the speed of light for very relativistic stellar models. For practical purposes, in these 
cases we just might as well set C s (r = R) = 0. Of course, this analysis does not hold for constant density models, but 
those are not considered in this paper, anyway. 

Let us now turn to the junction conditions at the surface of the star. We will always assume that e and C s go to 
zero when approaching the stellar surface r = R. Following the line of reasoning of we find that S is at least 
C 2 , whereas the differentiability of T depends on the value of p at the surface. For, if we let the subscripts in and 
ex represent the values for the interior and the exterior, respectively, we deduce from Eq. ( pl| ) that for the second 
derivative of T across the surface, the following relation has to hold: 



2A 



P\r=R 



(70) 



From condition (j^), it follows that the values of p at the surface depends on the value of e' at this point. As already 
discussed above, for polytropic equations of state, we have p = e = p' = at the surface. However, e' does not 
necessarily vanish, instead we have the following relation 



(71) 



which shows that the behavior of e' critically depends on the value of the polytropic index T. We can distinguish 
three different cases. For T < 2, we have e' — * 0, for T = 2, we have e' — > const., whereas for T > 2, we have e' — » — oo! 

This is somewhat disturbing, since for the boundary condition (^) this would mean that \p\ — > oo, unless the 
expression in brackets vanishes. However, this is not automatically guaranteed! Interestingly, the boundary condition 
( |69| ) for H is harmless for all values of T, since v' is always bounded. But, of course, in Eq. (^l|), we have to use p, 
which we have to compute from H using 



(P- 



J2-T 



H 



kT 



-H . 



(72) 



Again, we obtain an infinite value when T > 2, unless H vanishes at the surface. However, as with Eq. (|6£|), the 
boundary condition (|6^) does not guarantee the vanishing of H , even if H is initially set to zero. 

The reason for this peculiar behavior of p is that in the F > 2 case, the Eulerian description breaks down at the 
surface, which is a direct consequence of e' becoming infinite. In this case, a Lagrangian description of the fluid 
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perturbations would be much more appropriate, since for any polytropic equation of state, the Lagrangian energy 
density perturbation always vanishes at the surface, because of the vanishing of the Lagrangian pressure perturbation. 

What then happens to the Eulerian density perturbation? By definition, the Eulerian density perturbation Se is 
the difference between the total energy density and the background density e at the same location r, whereas the 
Lagrangian density perturbation Ae measures the density change in a fluid element that is displaced by some amount 
£*. It is through Taylor expansion to linear order that we obtain the connection between the Lagrangian and Eulerian 
perturbations: 

Ae = 8e + fV ■ (73) 

But of course, for the T > 2 case this expression is ill-defined at the surface, since the second term diverges. It is then 
clear that the Eulerian perturbation Se has to become infinite, too, in order to compensate for the blow up of e' and 
to yield a vanishing Lagrangian perturbation Ae. 

The whole discussion seems somewhat irrelevant, since we have seen that the equations can be written in terms of 
metric quantities only. However, we have to compute the second derivative of T, which, as can be seen fromEq. ©, 
depends on p. In the L > 2 case, we therefore must have a blow up of T" at the surface. Of course, this is very 
troublesome for the numerical discretization, and even for T = 2, we still have a discontinuity in T" , which can spoil 
the second order convergence of the numerical discretization scheme. 

The numerical evolutions indeed confirm the above analysis. By computing p with the aid of the Hamiltonian 
constraint (|57|), we find that for polytropic stellar models with L > 2, p tends to blow up at the stellar surface. For 
L < 2, we have p = at the surface, whereas for V = 2, we obtain a finite value p2]. 



IV. CONSTRUCTING INITIAL DATA 

To obtain physically valid initial data, we have to solve the constraint equations (^) and (|59|). As is well known, 
there is no unique way to do so, for those equations are underdetermined, and therefore one always has to make some 
additional assumptions about the geometry of the initial slice. 

The simplest way is to choose time symmetric initial data, since in this case the momentum constraints (^9|) can 
be satisfied trivially by setting all extrinsic curvature variables to zero, and we only have to solve the Hamiltonian 
constraint (p7[). Of course, time symmetric data are somewhat unphysical, since they represent a stage of a physical 
system which had an arbitrary amount of incoming radiation in the past. Nevertheless, they can give valuable insight 
into the behavior of the system under consideration. 

In general, all conceivable initial data for neutron star oscillations fall into two basic categories, or combinations 
thereof. First, we can assume the star being initially unperturbed. In the exterior vacuum region, we then can 
arbitrarily prescribe some metric perturbations representing a gravitational wave, which travels towards the star and 
eventually excites it to oscillations. The most convenient way to do so is to prescribe the initial Zerilli function Z 
and its time derivative Z, and to compute S, T, K, and K$ through Eqs. (|6^). In this way, we have initial data that 
automatically satisfy the constraints without having to solve any differential equation. 

The other possibility is to prescribe some fluid perturbations inside the star and then use the constraints to solve for 
the associated metric perturbations. As is well known, the matter distribution does not uniquely fix the metric, since 
we can always superpose some additional gravitational waves. It is therefore not clear at all what would represent 
"true" physically realistic initial data. 

A quite common procedure is to construct conformally flat initial data p4| . In our case, this amounts to setting 
S(t = 0) = So = in the Hamiltonian constraint j5^ ) and solving for To for a given initial fluid perturbation pq. 
Interestingly, the thus obtained initial data do not excite any of the spacetime modes to a significant extent, whereas 
initial data with the same fluid perturbation po, but with To = and So ^ 0, can show large w-mode excitations. 
In Fig. 1, we compare the wave signals of these two cases for time symmetric initial data. In the case where we set 
Sq = 0, there is no w-mode signal at all, and the waveform consists of pure fluid oscillations. In the other case, we 
can see a burst of gravitational radiation, which quickly damps away through the first w-mode. The final fluid ringing 
then coincides with the So — case. 

This shows that the presence of w-modes is not a generic feature of neutron star oscillations. They are only excited 
by a special subset of initial data. The open question now is whether or not they will show up in a real wave signal. 
Andersson and Kokkotas |^5[ have shown that by extracting the frequencies and damping times of the /-mode and 
the first w-mode in a gravitational wave signal, one can obtain important information such as mass and radius of 
the neutron star. These data, in turn, can then be used to restrict the possible equations of state. Of course, this 
method stands and falls with the presence of the w-modes in the wave signal. But as we have seen, it could well be 
that w- modes do not play a significant role at all. In a subsequent paper |p6|, we will investigate the excitations of 
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neutron star oscillations by means of particle scattering, which will show that the w-modes can only be excited by 
extremely relativistic particles. 

Let us now focus on initial data which are not time symmetric. In this case, we also have to solve the momentum 
constraints. Here, too, we have some additional freedom in specifying the extrinsic curvature. For instance we can 
choose TrK — 0, which is equivalent to setting 

K x = 4K 2 - 2K 5 , (74) 

or, equivalently because of Eq. (|f8|). 

r 2 K = 2r—l + 2(rX' + 2)K 2 - 3K 5 . (75) 
Or 

We then can substitute K in Eq. (|59| ) and obtain thus two coupled differential equations for the remaining unknowns 
K 2 and K§. For given fluid perturbations u\ and u 2 , the momentum constraints yield a unique solution. 

Before presenting the numerical results, we would like to comment on the applicability of the Lichnerowicz-York 
approach [ p7| for the construction of initial data for neutron star perturbations, which was used in [ pTfij ] . Stated in a 
somewhat crude way, this method consists in decomposing the six degrees of freedom of the extrinsic curvature into its 
trace TrK (1 DOF), a transverse trace free part A FTF (2 DOFs), and a trace free longitudinal part A^ TF (3 DOFs). 
Since the momentum constraints consist only of three equations with three source functions ji, which describe the 
matter distribution, only 3 DOFs are fixed, whereas the remaining ones are freely specifiable. The simplest way is to 
set TrK = Af FF = 0, i.e. the extrinsic curvature is purely longitudinal. It can thus be written in terms of a vector 
potential with three independent components, for which the momentum constraints provide a unique solution once 
the sources ji and appropriate boundary conditions have been specified. 

If the extrinsic curvature is first decomposed into spherical harmonics, this procedure still works. For the polar 
part, we then have four components and two constraint equations, i.e. only two components can be freely specified. 
Here, too, by assuming TrK — A FTF = 0, we can write the remaining longitudinal part in terms of a vector potential 
with two polar components. Solving the constraint equations with non- vanishing sources ji will yield nonzero values 
for those component, and therefore all four polar components of the extrinsic curvature will in general have nonzero 
values, either. 

But this is contradiction with the Regge- Wheeler gauge, which requires the component K4 to vanish. This could 
only be accomplished by setting one component of the vector potential to zero. However, this is not possible, since 
both components are already fixed by the constraint equations, and in general non-vanishing. 

Hence, if we want our initial data to be conform to the Regge- Wheeler gauge, we have to have K4 — 0, and therefore 
we cannot have both TrK — and Aj J TF = 0. By choosing TrK ~ 0, we already have used up all degrees of freedom 

to fix the extrinsic curvature components, since we have to set K4 = 0, too. 



V. NUMERICAL RESULTS FOR POLYTROPIC EQUATIONS OF STATE 



A. Some general remarks 

Our actual goal was, by using the ADM-formalism, to obtain the evolution equations in a hyperbolic first order 
form. However, due to numerical problems at the origin, we were forced to modify those equations in such a way that 
the resulting set (^) is more or less equivalent to the set of wave equations for the variables S, T and H. Thus, from 
the computational point of view, there is no point in sticking to the first order system, since the wave equations can 
be numerically integrated in a much faster way. Also, we will not explicitly integrate the fluid equations, since we can 
eliminate p in Eq. (|5l| ) by means of the Hamiltonian constraint (|57|) . 

For the numerical evolution, we therefore use Eqs. (|50| ) and ( |58D in the interior, whereas in the exterior we use 
Eqs. ((5^) and (|5l]) with p set to zero. Of course, in the exterior we could try to switch from S and T to the Zerilli 
function Z . We thus would have to integrate only one wave equation, and from the resulting waveforms, we could 
easily read off the emitted gravitational energy, since the radiation power is proportional to the square of Z. 

Unfortunately, this causes numerical problems, because at the seam, we would have to compute Z and its derivatives 
from S and T and vice versa. However, in order to obtain a correct value of Z , it is crucial that S and T satisfy 
the Hamiltonian constraint. On the other hand S and T satisfy the Hamiltonian constraint by construction, at least 
up to discretization errors, when directly computed from Z through Eqs. (|63|). At the seam somewhere outside the 
star, where we switch from S and T to Z ', we will always have some violation of the Hamiltonian constraint, thus the 
computed Zerilli function Z will not be quite correct. If we then go back and compute S and T from the just obtained 
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value of Z, they will differ significantly from their original values, since now they suddenly satisfy the Hamiltonian 
constraint. This mismatching at the seam gives rise to additional reflections, which rapidly amplify inside the star 
and cause the numerical code to crash after a few dynamical time scales. 

There is another point we would like to mention. When evolved with the Zerilli equation, the amplitude of an 
outgoing wave remains constant for large r. The same is true for S, however, the amplitude of T grows linearly with 
r. When computing Z from S and T by means of formula (0), this growth has to cancel, what is indeed the case 
when S and T exactly satisfy the Hamiltonian constraint. If not, Z will also start to grow for large r. But this is 
what happens, since numerically S and T do not satisfy the Hamiltonian constraint, ft is predominantly the high 
frequency components that get amplified the most strongly. Unfortunately, this can lead to quite rough waveforms 
for Z, while those for S and T look perfectly smooth. By increasing the resolution, this amplification will decrease, 
but the resolution has to be quite high in order to obtain accurate results. 

The most practical way to obtain a quite reliable Zerilli function, therefore, is not to completely switch to the Zerilli 
function in the exterior region, but to additionally evolve Z together with S and T. That is, close to the surface of the 
star we construct Z from S and T by means of formula ( pOj) and then use the Zerilli equation ( |6^ ) to independently 
evolve Z parallel to S and T. Of course, this amounts to the additional computational expenditure of evolving an 
extra wave equation, but we get rewarded by obtaining much more accurate results. 

To discretize the wave equations, we use the explicit second order leap-frog scheme. If the equations are correctly 
implemented, the numerical violation of the constraints should converge to zero in second order. By monitoring the 
Hamiltonian constraint in the exterior, this is indeed what we find. Of course, as mentioned above, for polytropes 
with T > 2, the stellar surface can reduce the convergence down to first order. 

We should make one final note concerning the numerical treatment of the origin r — 0. It is well known that in 
radial coordinates, the equations usually show a singular behavior close to the origin. Moreover, Taylor expansion 
around the origin shows that the equations admit two kinds of solutions, a regular and a divergent one. On physical 
grounds, one usually rejects the divergent solution. For the evolution, it is crucial that the numerical scheme preserves 
the regularity condition and suppresses the divergent solution. 

We have chosen the dynamical variables such that the only singular terms are the ones proportional to 1(1 + 1)/V 2 . 
The regularity condition requires all the perturbation variables to vanish at the origin, however, even for I — 2 the 
equations cannot be numerically evolved with a time step size At close to the maximal allowed time step size At max , 
which follows from the CFL-condition 

Ar 

At max = — , (76) 
c 

where Ar is the grid spacing in the r-coordinate and c the largest propagation speed, because the presence of the 
divergent 1(1 + l)/r 2 -terms causes the scheme to become unstable at the origin. Only by decreasing At can stability 
be reinforced. However, for large values of I, At has to be so small that it prevents one from obtaining numerical 
results within a reasonable time frame. 

We therefore propose a different way that allows stable evolutions with At close to At max for all values of I. Instead 
of having the boundary conditions located exactly at r = 0, we move it n grid points to the right, i.e. at r n = nAr, 
we impose S(r n ) = T(r n ) = 0. The actual value of n depends on I in the following simple way: n = I — 1. 

The numerical experiments and a detailed stability analysis |f22[ have shown that with this little trick, we can indeed 
obtain stable and second order convergent evolutions for arbitrary values of I with a Courant number of about 0.9. 



B. Results 

Results concerning the evolution of various initial data have already been presented in JllJ . Since the authors only 
focused on a single polytropic stellar model with a central density of eo = 3 • 10 15 g/cm 3 , we would like to consider 
three more models, one being less relativistic and the others being more relativistic. The physical parameters of the 
models are given in table 1. All results presented are for I = 2. Results for more stellar models and for other values 
of I can be found in |22| . 

As initial data, we choose a narrow time symmetric gravitational wave pulse, centered at ro = 80 km, where we 
prescribe T to have a Gaussian shape, and use the Hamiltonian constraint to compute S. The numerical resolution 
is 500 grid points inside the star. The resulting waveforms for the different models are shown in Fig. 2. Using a 
logarithmic scale, we plot the modulus of the Zerilli function, extracted at r — 100 km. 

The three waveforms clearly have quite different features. For the least relativistic Model 1 (upper panel), there is 
no w-mode signal at all, and the fall-off, which immediately follows the reflected wave pulse at about t = 0.6 ms, shows 
more a tail-like behavior. At t = 1.1ms it merges into the fluid ringing, which is dominated by the /-mode. This is 



14 



a somewhat unexpected result, since a direct mode calculation reveals that for this model there do exist w-modes, 
which should show up in the waveform. However, even the first w-mode has a quite large imaginary part (u> — 0.305 
+ 0.24i), which results from the model being less compact. Thus, the w-modes get buried in the tail- like fall-off. The 
more compact and therefore more relativistic the stellar model is, the less damped are the iw-modes, and they should 
eventually dominate over the tail-like fall-off. 

This is the case for Model 2 (middle panel), which is quite close to the stability limit. Here, instead of the tail-like 
fall-off, we now can see the ring-down of the first w-mode (oj — 0.256 + 0.085i). Still, it is much more strongly damped 
than the fluid modes, which therefore dominate the late time part. 

Model 3 (lower panel), which is unstable with respect to radial oscillations, is compact enough to allow the existence 
of trapped modes. Because of their comparatively weak damping, we cannot have a clear cut discrimination between 
the two ring-down phases, i.e. at first the w-modes and then the fluid modes. Instead, the signal consists of a mixture 
of both the trapped and the fluid modes. Nevertheless, the trapped modes dominate in early times, which can be 
seen from the Fourier transform in the upper panel of Fig. 3. In the lower panel, we have taken the Fourier transform 
at much later times, where the most strongly damped trapped modes have damped away and the fluid modes prevail. 
Still, the first two trapped modes are present in the spectrum. 

However, there is still one feature present in the wave signal for Model 3, which was also present in previous 
evolutions of the axial modes for very compact stars, but which apparently has not been noticed before (see e.g. Fig. 4 
of || and Fig. 5 of §). 

If, in the wave signal, we determine the frequency and damping time of the very first ring-down (from t ks 0.5 ms to 
about t k 0.9 ms), we find that they do not match to any of the quasi-normal modes associated with this particular 
stellar model. Instead, they are almost identical with the complex frequency of the first quasi-normal mode of an 
equal mass black hole! How is that? 

For black hole spacetimes, the quasinormal modes are derived from the Regge- Wheeler potential in the axial case 
and the Zerilli potential in the polar case. Both potentials exhibit a maximum at about r = 3M. For ordinary neutron 
stars the surface lies usually at a radius larger than 3M, therefore these peaks do not exist for those stars. However, 
for ultra-compact stars these peaks can lie outside. Any axial or polar incoming gravitational wave packet has to 
cross the corresponding peak before penetrating the neutron star. While trying to cross this peak, a large fraction of 
the wave will be scattered off, which in the black hole spacetime gives rise to the quasinormal ring-down, but since 
the situation is now almost identical for the neutron star spacetime, this associated ring-down will have black-hole 
characteristics. 

Since inside the star, the local speed of light e l/_M is largely reduced, it takes a while for the remaining wave packet 
to be reflected, and this is gives rise to the strong increase in the signal at about t = 0.9 ms. Here, the trapped 
wave packet finally finds its way out again, and it is only from thereon, that the wave signal consists of the proper 
fluid and spacetime modes. (It is also from this time when we have taken the Fourier transformation to obtain the 
upper spectrum in Fig. 3. If we had included the first ring-down, this spectrum would have sat on top of a large 
broad peak.) For less compact stars the local light speed inside the star is higher, and therefore the wave packet 
gets reflected much earlier, making the black-hole-like ring-down phase much shorter. We should mention that this 
ring-down phenomenon for very compact stars is much more pronounced for larger values of I p2[ . 

Of course, we gave a somewhat crude and heuristic explanation of this ring-down phenomenon; we therefore postpone 
a more detailed study to a future paper p9[| . 



TABLE I. List of the used polytropic stellar models and their physical parameters. 



Polytropic stellar models (F = 2, k = 100 km 2 ) 


Model 


eo [g/cm 3 ] 


M [M Q ] 


R [km] 


M/R 


1 


1.0-10 15 


0.802 


10.81 


0.109 


2 


5.0-10 15 


1.348 


7.787 


0.256 


3 


5.0-10 16 


1.031 


4.992 


0.305 
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VI. GETTING INTO TROUBLE: USING REALISTIC EQUATIONS OF STATE 



So far we have used polytropic equations of state, which are quite decent approximations to realistic equations of 
state as far as the bulk features of neutron stars like mass and radius are concerned. However, it is in particular 
the fluid oscillations that are very sensitive to local changes in the equation of state, which are due to the different 
behavior of the neutron star matter under varying pressure. It is therefore important to use realistic equations of 
state that take into account the underlying microphysics which determines the state of the matter as a function of 
pressure and temperature. 

Realistic equations of state cannot be given in analytic terms over the whole pressure range inside the neutron 
star, hence they usually exist in tabulated form only. To solve the TOV equations in this case, one has to interpolate 
between the given values in order to obtain the stellar model with continuous functions of radius r. In the following, 
we will make use of an equation of state called MPA |17|, which yields a maximal mass model of 1.56Mq. The 
non-radial oscillations modes for various stellar models have been compiled in fjCfl . 

If we try to repeat the evolution of the same initial data, but now for a realistic stellar model with, say, N — 200 
grid points inside the star, we will find, after a few oscillations, an exponentially growing mode that immediately 
swamps the whole evolution. It is clear that this has to be a numerical instability, since the frequency and e-folding 
time of this growing mode strongly depends on the chosen resolution. In fact, by further increasing the resolution, 
the growing mode starts to weaken and eventually vanishes completely. This happens at about N — 500 grid points 
inside the star, depending on the stellar model under consideration. Hence, we have this strange fact that for low 
resolutions the numerical evolution tends to be unstable, whereas for resolutions high enough, the evolution is stable. 
We should stress that this happens only with the use of realistic equations of state. It does not happen for polytropic 
stellar models, even if the resolution is extremely low. 

Further investigations show that the origin of this instability comes from the region close to the surface, where the 
sound speed has a sharp drop. In Fig. 4 we show the profile of the sound speed C s = y/dp/de for a stellar model 
using EOS MPA with a central density of eo = 4-10 15 g/cm 3 . It can be clearly seen that at r w 8.06km, there is 
a local minimum of the sound speed, where it drops down to C s rj 0.02. For larger r we can see a series of much 
smaller dips, which is an artefact of the numerical spline interpolation between the tabulated points. But the dip at 
r ~ 8.06 km is physical and is related to the neutron drip point, where the equation of state suddenly becomes very 
soft. Since this occurs still in the low pressure regime, the dip is present for any realistic equation of state. 

It is indeed this dip in the sound speed that is responsible for the numerical instability, for if we remove it "by hand" , 
we obtain a stable evolution! Moreover, the occurrence of the instability is independent of the actual formulation of 
the equations, since the terms that are responsible are the same in each case. But what are the "bad" ter ms ? B y 
simply crossing out individual terms, we find that the culprits are those terms in the fluid equations ( 52a ), ( 54a ), 
(|56|) or (58) which are not multiplied by C^. These are the terms which remain when the sound speed goes to zero, 
i.e. these which constitute the boundary conditions ( p8| ) and (|69|) . Without these terms, the evolution would be stable 
for any resolution. 

The next question is then whether this instability is just a feature of our finite differencing scheme and does not 
appear for a different scheme. It turned out that the instability does indeed occur for all discretization schemes we 
have tried: The central difference scheme for the wave equations, a central difference scheme on staggered grids and 
the two-step Lax-Wendroff scheme both for the first order system of equations. Of course, these are only explicit 
schemes, and we do not know whether switching to implicit schemes could be a remedy. 

However, in all above cases, the instability eventually vanishes when the resolution is increased above a certain 
threshold. For the central difference schemes, this threshold is at about N — 500 grid points, depending on the stellar 
model, whereas for the Lax-Wendroff scheme, the required minimum resolution in order to overcome the instability 
can even exceed N — 10000 grid points inside the star, which is clearly unacceptable. 

A more elaborate analysis of the nature of the instability can be found in J22|. There it is shown on the basis 
of a simplified toy equation that the resolution needed in order to overcome the instability strongly depends on the 
location and depth of the dip in the sound speed. The closer to the surface and the deeper the dip, the higher the 
required resolution. Since this resolution is only needed in a very tiny region close to the surface, it would be enough 
to just refine the grid in this region and use a coarser grid outside. This could be accomplished, for instance, by fixed 
mesh refinement, since this region is determined by the profile of the sound speed, which does not change throughout 
the evolution. However, we then would have to deal with the transition from the coarse grid to the fine grid and 
vice versa, which might be troublesome. Another drawback is that for each different stellar model, we would need a 
different grid refinement, and it would be a matter of trial and error to find the appropriate refinements for a stable 
evolution. 

Yet, there is a better way out. We can try to find a new radial coordinate x, which is related to the actual radial 
coordinate r in such a way that an equidistant grid in x would correspond to a grid in r that becomes automatically 
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denser in regions where the sound speed assumes small values. A simple relation between the grid spacings Ax and 
Ar that has the desired property is 

Ar(r) = C s (r)Ax . (77) 

An equidistant discretization with a constant grid spacing Ax would correspond to a coarse grid in r for large values 
of C s , which becomes finer and finer as the sound speed C s decreases. 

From the relation (f77|), we can immediately deduce the form of our new coordinate x as a function of r. By replacing 
the A's by differentials we obtain 

dr C s (r) 

or 

x{r) = .[<m- (79) 

As a consequence, the derivatives transform as 
and 

d 2 2 d 2 , d 

— , = C s —7T + C S C S — . (81) 
dx 2 dr 2 dr 

From the last relation, we see that the thus defined coordinate transformation will transform the wave equation in 
such a way that the propagation speed with respect to the x-coordinate will be unity throughout the whole stellar 
interior. 

Of course, we have to use relation (JT^) with some caution, for if C s — this transformation becomes singular. 
And this is what happens at the stellar surface. If, for instance, the profile of the sound speed is given in the form 
C s — C*o(l — r/R), where we have C S (R) = 0, we can find an analytic expression for x 

s(r) = -^|]og(l-~) , (82) 

which tells us that at the surface r — R we obtain x(R) = oo. In this case, the coordinate transformation seems to 
be quite useless, since numerically we cannot deal with a grid that extends up to infinity. We would have to truncate 
it somewhere. But from a numerical point of view, this is not that bad, since going to infinity in the x-coordinate 
would correspond to an infinitely fine resolution in the r-coordinate at the stellar surface. But this is numerically 
impossible as well, so truncating the x-coordinate at some point means to define a maximal resolution in r at the 
surface. However, this case usually does not happen. In all cases considered, x(R) has always a finite value. 

It should be noted that the above transformation is of the same kind as the definition of the tortoise coordinate r* 
in Eq. (^), which also leads to wave equations with unity propagation speed throughout the whole domain. Thus, 
we may call x a hydrodynamical tortoise coordinate. 

Figure 5 shows the r-coordinate as a function of the x-coordinate for three different stellar models. From the curves, 
it is clear that an equidistant grid spacing in x will result in an equivalent spacing in r that gets very dense towards 
the surface of the star, which means that this part gets highly resolved. And this is exactly what we need in order to 
overcome the instability and to obtain a decent accuracy. 

To obtain the background data on the x-grid, we have to solve the TOV equations (||) with respect to x. Since we 
also need r as a function of x, we simultaneously solve Eq. (uaj), too. The transformed set of equations then reads: 



£ = -r<p+-o < 83<! > 

ax ax 

(IV 

^ = C S . (83d) 
dx 
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Of course, this transformation is only defined in the stellar interior for the fluid equations, for it is only there that 
the propagation speed is the speed of sound (apart from the factor e"~ A ). If we still wanted to use Eq. (58), where in 



the interior T plays the role of the fluid, we would have to switch from the x-grid in the interior to the r-grid in the 
exterior, which is somewhat inconvenient. It is therefore much more natural to explicitly include the fluid equation 
(|56|), which is defined in the interior only. Thence, it is only Eq. (|5^) which will be transformed according to the 
transformation ((78]), whereas we keep the wave equations for S and T as they are given in ( |50|) and (|5l|). 
The transformed fluid part of the fluid equation <Mf) reads 



d 2 H 



d 2 H ( , , . v a C sx \ dH 

+ (2v x + \ tX ) 



dx 2 V C's Cs J dx 



r 



rC K rC, 



(84) 



Here, the subscript x denotes a derivative with respect to x. The last missing thing is the transformation of the 
boundary condition (|6|). Unfortunately, we cannot transform Eq. ( p9| ) in a straightforward way, since at the surface 
it is C s = 0, and therefore the transformation of the derivative d/dr — C~ 1 d/dx is not defined. However, the inverse 
transformation ( |80| ) can make sense if we note that if C s = then any derivative with respect to x has to vanish. This 
is in particular true for H itself. Thus, at the stellar surface r = R, where the sound speed C s vanishes we impose 
the following boundary condition for H: 



dH 

dx 



X(R ) = . (85) 



This corresponds to reflection at a loose end on the x-grid. 

Since we only transform the fluid equation and not the equations for the metric perturbations S and T, we have 
to simultaneously use two different grids: The x-grid in the interior for the fluid variable H, and the r-grid both in 
the interior and exterior for the metric variables S and T. Because of the coupling, at each time step we have to 
interpolate H from the x-grid onto the r-grid in order to update equation (plf), and, vice versa, both S and T from the 
r-grid onto the x-grid in order to update equation fl5q) . This can easily be accomplished by using spline interpolation. 

In Fig. 6 we demonstrate the convergence of this new method using a stellar model based on the realistic EOS MPA. 
For initial data similar to the ones used for Fig. 2, we plot the evolution of a Gaussian wave packet together with 
the evaluation of the Hamiltonian constraint for three different spatial resolutions. By doubling the resolution, the 
error in the Hamiltonian constraint should decrease by a factor of 4, if the scheme is of second order. Conversely, the 
magnitude of the error should remain roughly constant, if for each doubling of the resolution, the error is multiplied 
by a factor of 4. This has been done in Fig. 6, and one can clearly see that the magnitudes of the thus rescaled errors 
remain constant for the different resolutions. Furthermore the scheme is stable since the error does not grow with 
time. 

The fluid modes have their largest amplitude in the region close to the surface. To obtain the correct mode 
frequencies in the power spectrum, it is important to have high resolution close to the surface of the star, and this 
is what is accomplished by our new coordinate x. In Fig. 7, we show the power spectra of waveforms obtained from 
runs with different resolutions. As initial data, we chose a narrow fluid deflection at the center of the star together 
with 5 = 0. In this case, we suppressed any w-mode contribution. 

It is apparent that even for the quite moderate resolution of 50 grid points on the x-grid inside the star we obtain 
very accurate frequencies for the first couple of fluid modes. On the r-grid, however, a resolution of 200 grid points 
is still not enough to obtain the same accuracy, and the peaks of the higher p-modes in the spectrum are still quite 
far off their true values. Of course, these results were obtained with a polytropic equation of state, since otherwise 
we would not have been able to perform the evolution on the r-grid at all, because of the occurring instability. 

Finally, we want to demonstrate the effectiveness of our coordinate transformation by evolving an initial perturbation 
for the MPA equation of state. We choose the central energy density to be 3-10 15 g/cm 3 , yielding a stellar mass of 
AI = 1.49M Q . The resolution is N = 200 grid points, which, again, would be too low to yield a stable evolution using 
the r-coordinate. Here we do not have any problems, the evolution is stable for any chosen resolution. In Fig. 8, 
we show the power spectrum of the evolution of a sharp initial peak in the fluid perturbation H , which leads to the 
excitation of a multitude of fluid modes. In the interval up to 100 kHz, we can find 37 modes. The direct calculation 
of the first couple of modes with an eigenvalue code shows perfect agreement, which is demonstrated in the lower 
panel of Fig. 8. Actually, the direct mode calculation is much more involved and time consuming than the evolution 
with a subsequent Fourier transformation. Even with a thousand grid points inside the star, it is a matter of minutes 
to obtain an accurate fluid mode spectrum for stellar models with realistic equations of state. Of course, with this 
method one can only obtain the real parts of the frequencies, the imaginary parts are too small to be determined. 



18 



VII. CONCLUSIONS 



In this paper, we have derived, using the ADM-formalism, the perturbation equations for non-rotating, spherically 
symmetric neutron star models with barotropic equations of state. We have shown how to choose shift and lapse and 
the initial data on the three dimensional initial slice in order to obtain the Regge- Wheeler gauge. We pointed out 
that not only some of the metric components but also one of the extrinsic curvature components has to vanish in 
order to preserve the Regge- Wheeler gauge throughout the evolution. 

Unfortunately, the equations that directly come out after the expansion in spherical tensor harmonics are not suited 
for the numerical treatment due to problems at the origin. We therefore had to introduce different variables that 
remove the singular structure of the "raw" equations. Since the final system of equations in those variables is nothing 
else but a system of coupled wave equations written in first order form, it is therefore much more convenient to 
use instead the second order wave equations for the numerical evolution. Those wave equations have already been 
previously derived by Allen et al. in a different way. Still, our formulation is quite useful because it facilitates the 
construction of initial data in terms of metric and extrinsic curvature perturbations. 

Moreover, by introducing the ADM-formalism for the derivation of the evolution equations of non-rotating neutron 
stars, we have performed the first step towards the much more involved task of deriving the evolution equations for 
rotating neutron stars. Here, the linearization of the original Einstein equations yields quite obscure perturbation 
equations, which first have to be cast into a suitable hyperbolic form before trying to discretize them on a numerical 
grid. With the ADM-equations this is a much more straightforward procedure. 

We have presented results for different stellar models with polytropic equations of state. For the least relativistic 
model, the wave signal did not show any w-modes, although mode calculations show that they do exist. But they 
are so strongly damped that the tail-like fall-off wins over the w-modes. Only for more compact stellar models, 
the w-modes are less strongly damped, and they can show up in the waveforms. However, even for very relativistic 
neutron star models, w-modes are not necessarily part of the waveform. In fact, their occurrence strongly depends on 
the choice of initial data. For initial data representing a gravitational wave pulse that will hit the star, w- modes in 
general will be part of the spectrum, however, for initial fluid perturbations, the choice of conformally flat initial data 
(So = 0) can totally suppress the presence of w-modes, whereas setting T = and Sq ^ 0, which can be regarded in 
some sense as a maximal deviation from conformal flatness, shows strong w-mode excitation. 

For very compact stars, we found that the complex frequency of the very first ring-down phase does not correspond 
to any of the proper quasi-normal modes of the stellar model. Instead, the frequency is almost identical with the first 
quasi-normal frequency of a black hole with the same mass. We explained this phenomenon by observing that those 
parts of the potential, which are responsible for the first scattering of the impinging wave packet, lie outside the star 
and therefore are the same for a black hole. Only if the part of the wave packet that could penetrate the neutron 
star is reflected, the proper ring-down with the trapped and fluid modes starts. This phenomenon, although already 
present in previous evolution results [p8p|l, apparently has not been recognized before. 

When switching to realistic equations of state, we are faced with a numerical instability that is due to the sharp 
drop of the sound speed at the neutron drip point, which is located quite close to the stellar surface. Since the 
neutron drip point is located in the low pressure regime, it is a feature of any realistic equation of state, and therefore 
the instability should occur for all of them. The quite peculiar fact about this instability is that it disappears when 
the resolution is increased above a certain threshold. This threshold depends on the depth of the dip of the sound 
speed and is the higher the deeper the dip is. It seems that this instability is present for any explicit numerical 
discretization scheme, however, the minimum resolution that is needed in order to suppress it, can vary drastically 
for different schemes. Indeed, when using the two step Lax-Wendroff scheme for the equations written as a first order 
system, we were confronted with an almost unsurmountable resolution threshold needed for stable evolution. 

This instability does not depend on the actual formulation of the equations and is present for any of the given 
formulations of this paper. In fact, it is not only a peculiar problem of the non-radial perturbation equations, the 
same instability also occurs in the radial case ]22] ] ; the reason being that the structure of the fluid equations is the 
same in both cases. 

As a remedy, we introduced the coordinate transformation (78), which leads to a wave equation for the fluid with 
unity propagation speed throughout the whole star. It is now possible to evolve oscillations of realistic stellar models 
for any resolution. As a further benefit, the evolution of the fluid on the x-grid is, for a given number of grid points 
inside the star, much more accurate than the corresponding evolution on the r-grid with the same number of grid 
points. This is because the ^-coordinate resolves the region close to the stellar surface much better than the ordinary 
r-coordinate. 

This coordinate transformation could also prove useful for fast rotating neutron stars, since the evolution equations 
of the fluid will certainly have the same structure that leads to the numerical instability in the non-rotating case. 
Since the rotating case is a 2D problem, resorting to higher resolutions may become impossible, and therefore our 
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coordinate transformation could do the job. Maybe, this transformation could even be of avail in nonlinear evolutions, 
although it is not that straightforward to implement, since there is no fixed background model with a static profile 
of the sound speed. A possibility is to view Eq. ( |78| ) as a dynamic coordinate condition that depends on the actual 
profile of the sound speed at every instance of time. 
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FIG. 1. Wave signals from two different sets of initial data. The signal coming from initial data with To = (solid line) 
shows a strong gravitational wave burst followed by a w-mode ring-down, whereas the signal coming from data with So = 
(dashed line) does not show any w-modes at all but rather consists of pure fluid ringing. Both signals coincide after the to-mode 
has damped away. 
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FIG. 2. Time evolution of the Zerilli function for the three different polytropic stellar models Ml, M2, and M3 with I = 2. 
The numerical resolution is 500 grid points inside the star. In the least relativistic model (Model 1, upper panel), there are 
no w-modes observable. For Model 2 (middle panel), both the w- and the fluid modes are excited. Model 3 (lower panel), 
which is the most compact, exhibits two ring-downs, one ranging from 0.6 to 0.9 ms and corresponding to the first black hole 
quasi-normal mode, and the "proper" ring-down starting at about 1.0 ms and consisting of the various trapped and fluid modes. 
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FIG. 3. Power spectra for Model 3 with 1 = 2. In the upper panel, the Fourier transform is taken for an early starting time 
(t > 0.9 ms) and shows the presence of the weakly damped trapped modes. In the lower panel, the Fourier transform is taken 
at a much later time, where most of the trapped modes have damped away and the p-modes prevail. Still, the first two trapped 
modes are clearly visible. 
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FIG. 4. Profile of the sound speed C s inside the neutron star model using EOS MPA, and a section near the surface (inlet). 
At the neutron drip point around r ^ 8.06 km, the sound speed exhibits a local minimum. 
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FIG. 5. The r-coordinate as a function of the ^-coordinate for three different stellar models. An equidistant discretization 
of x corresponds to a increasingly finer resolution in r as the stellar surface is approached. 
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FIG. 6. It is shown, using a realistic EOS, the evolution of the metric variable T together with the error in the Hamiltonian 
constraint (HC) for the three different resolutions of N = 200, N = 400 and N = 800 grid points inside the star. The graphs 
for T practically coincide for the different resolutions. In doubling the resolution, the error should decrease by a factor of 4, 
which is compensated by multiplying the error by a factor of 4 for N = 400 and by a factor of 16 for N — 800. It is clearly 
discernible that the magnitudes of the scaled errors remain constant, which shows the second order convergence of our new 
numerical scheme. In addition the errors do not grow with time, which shows the stability. 
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FIG. 7. Comparison of the power spectrum obtained from an evolution with a resolution of N x = 50 grid points on the 
rc-grid with the spectra obtained from evolutions with different resolutions on the r-grid. Even for only 50 grid points in x, the 
peaks are closer to the true values than for 200 grid points in r. 
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FIG. 8. Upper panel: Power spectrum of the evolution of a sharp initial fluid pulse with So = on the x-grid. No w-modes 
are present, but dozens of p-modes are excited. In the lower panel, we include the modes computed by explicit mode calculation. 
Note that in contrast to polytropic equations of state, the modes do not lie at equidistant locations. 
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